El Nino and the Delayed Action Oscillator 
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We study the dynamics of the El Nino phenomenon using the mathematical model of delayed- 
action oscillator (DAO). Topics such as the influence of the annual cycle, global warming, stochastic 
influences due to weather conditions and even off-equatorial heat-sinks can all be discussed using only 
modest analytical and numerical resources. Thus the DAO allows for a pedagogical introduction 
to the science of El Nino and La Nina while at the same time avoiding the need for large-scale 
computing resources normally associated with much more sophisticated coupled atmosphere-ocean 
general circulation models. It is an approach which is ideally suited for student projects both at 
high school and undergraduate level. 
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I. INTRODUCTION 

Where weather is concerned, attempts at its predic- 
tion [1, 2] and its apparent unpredictability have al- 
ways intrigued the human mind. Related phenomena 
such as storms, floods and droughts, as well as global 
weather and issues of climate prediction such as the 
green-house-gas effect and global warming regularly ap- 
pear in news bulletins and TV programs worldwide [3, 4]. 
Thus weather is where many of us get our first expo- 
sure to science and in particular science's attempts at 
predicting natural behavior of highly complex systems. 
However, precisely because of this complexity, inexperi- 
enced science enthusiasts rarely have the tools needed to 
engage in predicting atmospheric behavior. 

One particular example of a natural atmospheric phe- 
nomenon which continues to attract regular public at- 
tention is the so-called El Nino event in the equato- 
rial Pacific, which occurs irregularly at intervals of 2-7 
years with occasionally dramatic consequences for many 
regions worldwide. These "teleconnections" can be ob- 
served in many places, for example spring rainfall levels 
in Central Europe, flooding in East Africa, and the fe- 
rocity of the hurricane season in the Gulf of Mexico [5] . 
Predictions of when and where El Nino weather patterns 
will dominate local weather conditions are very hard to 
obtain. Usually, so-called coupled atmosphere-ocean gen- 
eral circulation models (CGCMs), combined with large 
scale computing resources are needed to study and pre- 
dict El Nino's consequences [6, 7]. Clearly, such resources 
will normally be beyond the reach of teachers and even 
lecturers at university. Consequently, student projects 
on the El Nino phenomenon which involve a substantial 
quantitative component are rare. 

Fortunately, a route towards qualitative and quan- 
titative prediction exists. This approach is based on 
the delayed-action oscillator (DAO), a one-dimensional, 
non-linear ordinary differential equation with delay term 
which was introduced in 1988 by Suarez and Schopf [8]. 
The DAO models the irregular fluctuations of the sea- 
surface temperature (SST), incorporating the full cou- 



pled Navier-Stokes dynamics of the El Nino event via a 
suitably chosen non-linearity. While previous studies of 
the model have concentrated on the universal properties 
of the dimensionless dynamics [8, 9], we shall show in the 
present paper how to use the DAO to predict many of the 
essential features of El Nino events. We shall extend the 
model to quantitatively include the effects of the annual 
cycle, global warming, randomness in temperature condi- 
tions and the possible influence of the dynamics outside 
the equator. We believe that this discussion of the DAO 
opens up various avenues for student projects, be it at 
the high school or undergraduate level. Important un- 
derlying physics and meteorology concepts such as ocean 
waves and their speeds, non-linear dynamics, the Corio- 
lis force, convective cycles of ocean and atmosphere cur- 
rents, global wind patterns and their stability can be dis- 
cussed and help the students understand that El Nino — 
far from only being a bringer of doom as sometimes por- 
trayed in the media — is in fact a natural consequence 
of our planet's celestial dynamics [10]. 

The paper is organized as follows. In Section II, we re- 
view the El Nino-Southern Oscillation phenomenon. Sec- 
tion III introduces the basic DAO model and we study 
its basic dynamics. Sections IV, V and VI investigate the 
influences of annual forcing, global warming and random 
effects on the DAO, with Section VII incorporating all 
these processes. Section VIII then extends that single 
DAO to more than one coupled DAO. We conclude in 
Section IX. Information on the numerical methods and 
an example code can be found in the appendix. 



II. ENSO PHENOMENOLOGY 

The term El Nino was originally used, from as early 
as the late 19th century [11], in reference to an observed 
warm southward flowing current that moderates the low 
SST along the west coast of Peru and Ecuador [12, 13]. 
This current arrives after Christmas, during the early 
months of the calendar year, and was therefore named 
"El Nino" (The "Little Boy" or Christ Child in Spanish). 
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In certain years this current would be unusually strong, 
bringing heavy rains and flooding inland, but also deci- 
mating fishing stocks, bird populations and other water- 
based wildlife in what would normally be an abundant 
part of the Pacific. 

Today, the term El Nino is most often used when de- 
scribing a far larger-scale warm event that can be ob- 
served across the whole of the Pacific Ocean by cer- 
tain characteristic climatic conditions. El Nino is just 
one phase of the so-called El Nino-Southern Oscillation 
(ENSO) phenomenon, i.e. an irregular cycle of coupled 
ocean temperature and atmospheric pressure oscillations 
across the whole equatorial Pacific. In Fig. 1, we show 
the measured SST as well as the average annual cycle. 
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FIG. 1: Observed ENSO SST oscillations since 1980. The 
average (dashed black, right axis) and measured (grey, left 
axis) SST's, minus 27.1°C, in the Nino 3.4 region (cp. Fig. 2) 
with the SST anomaly (solid black, left axis). El Nino events 
are denoted by a black square above the event, La Nina events 
(see section II C) by a black circle below. 

Since 1899, there have been 29 recorded El Nino events, 
giving an average period for ENSO of 3.7 years [14] . How- 
ever, there is evidence to suggest that El Nino events are 
becoming stronger and more frequent; there have been 6 
events in the last 20 yrs (an average period of 3.3 years), 
2 of which were the strongest of the century [12, 14]. 
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FIG. 2: Schematic view of the Nino regions 3 (right hashed 
rectangle), 4 (left hashed rectangle) and 3.4 (rectangle with 
dashed line borders). The thermocline height in normal 
(solid) and El Nino conditions (dashed line) is indicated 
schematically in the diagram. 



The SST is about 8°C higher in the West, with cool tem- 
peratures off South America, due to cold water upwclling, 
i.e. rising from the deeper layer to the surface [15]. This 
cold nutrient-rich water supports diverse marine ecosys- 
tems and major fisheries, so upwelling is of immense eco- 
logical and economic importance to the west coast of the 
Americas [12, 13]. 

The warmer SST in the West means increased evapora- 
tion and rising moist air, and hence high levels of rainfall 
(a tropical climate), along with low atmospheric sea-level 
pressures in the West. In the East, the colder SST re- 
sults in far less rainfall, an arid inland climate and higher 
atmospheric sea-level pressures. This induced pressure 
difference drives the circulation of the trade winds. A 
measure of the strength of the trade winds is the South- 
ern Oscillation Index (SOI), defined as the normalized 
difference in sea-level pressures between Tahiti and Dar- 
win, Australia. 

The whole process is a therefore a self-sustaining cycle; 
winds cause thermocline and SST changes, which cause 
pressure differences, which in turn drive the winds. 



Normal conditions 



El Nino conditions 



In order to understand the El Nino anomaly, let us 
first describe the normal conditions that exist in the Pa- 
cific Ocean between 120° East and 90° West. In Fig. 2 
we show a schematic view of the relevant section of the 
Pacific. 

In normal conditions, the trade winds blow from East 
to West across the tropical Pacific, resulting in warm sur- 
face water being pushed west, so that the thermocline — 
a region of water roughly 100m below the surface where 
there is a very sharp temperature gradient separating the 
warm surface layer from the cold abyss — is depressed 
further in the West, and made shallower in the East [11]. 



During an El Nino event the SST and SOI change dra- 
matically. In the central and western Pacific, there is 
relaxation in the strength of the trade winds. This leads 
to a depression of the thermocline in the East and an 
elevation in the West, because less warm water is be- 
ing pushed west. This increase in thermocline depth in 
the East can be as large as a factor of three, leading to 
an increase in SST, and the effective cut-off of upwclling 
of nutrients to the surface. This has devastating conse- 
quences for the primary producers (such as plankton) in 
the food chain, with knock-on effects for other life, such 
as fish and birdlife [11]. 
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The increased SST in the East also leads to an in- 
creased evaporation, and large amounts of rain, often 
causing extensive flooding. The rising air decreases the 
normally high pressure. During large El Nino events the 
SOI is very negative, indicating a reversal in trade winds. 
This causes a great deal of warm, moist air (and hence 
rain) to move east, whilst leaving a severe drought in the 
western Pacific, often resulting in large bush fires. 

C. La Nina conditions 

La Nina (The "Little Girl") is the phase of the ENSO 
cycle which is opposite to El Nino. It is characterized 
by unusually cold ocean temperatures in the central and 
eastern Pacific. La Nina arises by a strengthening in the 
trade winds, further depressing the thermocline in the 
West (and elevating it in the East), which increases cold 
water upwelling in the eastern Pacific. This lower SST 
means the high levels of nutrients can support an abun- 
dance of life in the Ocean. The SOI is also positive (in- 
creased pressure difference), bringing enhanced rainfall 
in the West. 



III. THE DELAYED-ACTION OSCILLATOR 

A. A non-linear oscillator with delay term 

The DAO is a simple non-linear model in which the 
ENSO effect is modeled by just three terms [8]. With 
T denoting the temperature anomaly, i.e. the deviation 
from a suitably defined long term average temperature, 
the model is given as a first-order, non-linear differential 
equation 

dT 

— = kT- bT 3 - AT(t - A) (1) 

with appropriate coupling constants k, b, A and A. The 
first term on the right hand side in (1) represents a strong 
positive feedback within a coupled ocean-atmosphere sys- 
tem which is assumed to occur in the Nino 3.4 region. A 
much deeper thermocline as in the West, or a cold SST 
as in the East, would result in much poorer feedback 
coupling. The feedback process works through advective 
processes giving rise to temperature perturbations which 
in turn result in atmospheric heating. This heating leads 
to surface winds, driving the ocean currents which will 
then enhance the anomaly T. 

The second component is an unspecified nonlinear net 
damping term that is present in order to limit the growth 
of unstable perturbations. This term represents various 
ocean advective processes and moist processes that stop 
temperatures diverging. 

Lastly, the model also considers equatorially-trapped 
ocean waves propagating across the Pacific, before inter- 
acting back with the central Pacific region after a cer- 
tain time delay. These ocean waves are "hidden" Rossby 



waves which move westward on the thermocline, reflect 
off the rigid continental boundary in the West and then 
return eastward along the equator as Kelvin waves. Their 
respective wave velocities and the width of the Pacific 
basin determine the transit time, which is denoted A in 
Eq. (1). The strength of the returning, emerging signals 
relative to that of the local non-delayed feedback is de- 
noted by A. We also assume that < A < k due to loss 
of information during transit and imperfect reflection. 
We shall often use the DAO in its dimcnsionless form 

dT' 

W = T'- (T') 3 - aT'(t'-S), (2) 

where we ha ve scaled time by t' = kt and temperature by 
T" = \JbjkT . The new constants then are a = A/k and 
S = kA. In what follows, we shall for simplicity suppress 
the prime in all expressions where there is no ambiguity. 

B. DAO wave propagation 

The DAO delay term has a negative coefficient rep- 
resenting a negative feedback. To see the reason for 
this, let us consider a warm SST perturbation in the 
coupled region. This produces a westerly wind response 
that deepens the thermocline locally (immediate positive 
feedback), but at the same time, the wind perturbations 
produce divergent westward propagating Rossby waves 
that return after time A to create upwelling and cool- 
ing, reducing the original perturbation. This process can 
therefore be seen as a phase-reversing reflection off the 
coupled region; the resultant upwelling will cause further 
Rossby wave propagation, which in turn will cause down- 
welling and warming after a further delay, A. Hence we 
expect the period of model oscillations to be no less that 
2 A — "period doubling" [8]. 

Note that the eastern boundary will also contribute a 
delay term to the DAO, however equatorial Kelvin wave 
reflection from eastern boundaries is weak as the bulk of 
the wave propagates along the coast [16]. For this reason 
we may ignore these effects. 

It is possible to estimate the wave transit times, and 
hence the time delay A of the DAO model, from knowl- 
edge of the wave speeds for Rossby and Kelvin waves 
(see Appendix A). We shall return to this direct com- 
parison later, but for now we focus on the behavior of 
the dimensionless DAO (2). 

C. Linear stability analysis 

Equation (2) cannot be solved analytically. But wc 
may study its dynamics by analyzing its behavior close 
the its fixed points, i.e. those special solutions T of (2) for 
which dT/dt = which also immediately implies T(t) = 
T(t — S) for all t. Hence we need to solve 

(1 -a)T- f 3 = (3) 
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and so the fixed point solutions are T = 0, ±\A — ol. 
To = is the trivial fixed point, corresponding to a 
zero temperature anomaly. It is also unstable, i.e. any 
small perturbation will lead to a system moving away 
from 0. Therefore we shall concentrate on the other two 
fixed points T± = ±\/l — a. Due to symmetry, it is fur- 
thermore sufficient to consider one of these only and we 
choose f + = VI — a - 

The dynamics close to T + can be examined by a linear 
stability analysis. Let us consider a small perturbation 
S = T — f + . Ignoring terms of more than linear order in 
S [8], we find when substituting S back into (2) that 
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7 Q 

— = (3a-2)S- aS(t-5) 



(4) 



We now look for solutions of the form S = Ae at , where 
A is a constant real number equal the initial condition 
S(0) and a is a complex number. Substituting this in (4) 
gives an equation for a 



a = (3a — 2) — ae 



-aS 



(5) 



Writing a = ctr + ia\, where ctr, cti are both real, we find 

<7 R = (3a - 2) - a cos(ct i( 5) e~ aR5 (6a) 
cri = a sin(cri<5) e- a * s (6b) 

When ctr = 0, the solution S = Ae lait is purely oscil- 
latory and hence the fixed point is a center — neither 
stable nor unstable — with a finite frequency of 2n/ai. 
These neutral solutions occur along (an infinite number 
of) neutral curves of the form: 



8 = 



2mr =F arccos (^-^) 
^a 2 -(3a-2f 



n e N U {0} (7) 



In Figure 3, we show the (a, 8) phase diagram of the 
DAO obtained from the linear stability analysis. Only 
0.5 < a < 1 is plotted, since for this domain the coeffi- 
cient on the local term in (4) is always smaller than that 
of the delayed term. If it were larger, then the local term 
would always dominate, and the solution would always 
be stable, i.e. for a < 0.5 the fixed point is always stable 
and attracting, regardless of the initial conditions. 

Choosing (a, 8) below the first neutral curve, e.g. 
within the shaded region of Figure 3, corresponds to 
ctr < and hence the solution is a decaying spiral - 
the fixed point is stable and attractive. However, tak- 
ing (a, 8) above this neutral curve gives ctr > 0, and so 
the perturbation spiral grows outwards — the fixed point 
is unstable and no matter how close to T the solution 
starts, it will always move away. 

If (a, 6) is taken above the second neutral curve, the 
fixed point is still unstable (ctr > 0) but the number of 
unstable modes has increased to two. However, the basic 
DAO equation (2) is only a one-dimensional problem, as 
there is only one variable T, and so the extra unstable 
modes do not influence the dynamics. 
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FIG. 3: Phase diagram for the stability of fixed points T+. 
Neutral curves for n = (solid line) and n = 1 (dashed line) 
are shown together with stable (grey) and unstable regions 
(white) . 



This linear stability approach is highly useful, given 
that it depends on only two parameters, a and 8, and 
reveals valuable information about the stabilities of the 
fixed points. Given an a and a 8, it also provides a 
solution to the DAO given by: 



T(t) =T++ Ae" Rt e 



%<7\t 



(8) 



which is perfectly valid to within an order of 0(T), pro- 
vided it stays close to T+. However, once a solution is too 
far from an unstable fixed point or if the initial condition 
starts too far away from the fixed point, this linear stabil- 
ity approach breaks down. Figure 4 shows linear stability 
solutions for different 8 values, given by Eq. (8), along 
with the full solutions discussed in Section HID. 



D. Computation of full solutions by numerical 
integration 

Let us now construct solutions to the DAO which are 
outside the region of validity of the linear stability anal- 
ysis. As mentioned before, this is no longer possible 
analytically, so we will use numerical techniques to in- 
tegrate up the time-behavior of (2). Fortunately, this 
involves the use of standard numerical tools available in 
a variety of mathematical software packages. For con- 
venience, we have chosen the predefined Mathematica 
package NDelayDSolve [17]. The Mathematica code used 
to solve the basic DAO equation (2) is part of the code 
shown in Appendix B. 

In Figure 4, we show the results obtained from Eq. (2) 
when carrying out this numerical procedure for a = 0.75 
and various values of 8. The initial value used is T(0) = 
0.55 which is close to the fixed point T + = 0.5. In Fig. 
5, we show the (T,dT/dt) phase portraits for two values 
of 8 close to the first neutral curve. 

Let us now briefly discuss the observed behavior. In 
Fig. 4(a), we see that the numerical solution and the 
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FIG. 4: An a = 0.75 cross-section for T vs t for (a) a stable 
solution with 5 = 1, (b) a nearly neutral solution at S = 2 
and (c) an unstable oscillation at 5 = 4. Solid lines indicate 
numerical solutions for the anomaly T, dashed lines are ob- 
tained from the linear stability solutions T+ + S. For clarity, 
we only plot the 5 = 4 linear stability solution up to t = 30. 



linear approximation arc in good agreement and quickly 
converge to the stable fixed point T + = 0.5 as expected 
for 5 = 1. In Fig. 4(b), we are still close to the neutral 
curve but have already crossed over into the unstable 
regime. Consequently, the numerical solution starts to 
slowly spiral away from T + = 0.5, becoming large enough 
so that the dynamics are also influenced by the second 
unstable fixed point at f_ = —0.5, resulting in the onset 
of oscillations between them. This is shown in more detail 
in Fig. 5(b). The agreement with the numerical solution 




FIG. 5: Phase portraits of vs T with a = 0.75, solved up 
to tfinai = 1005 with (a) <5 = f.6 and (b) 5 = 2. Solid lines in- 
dicate numerical solutions of Eq. (2), dashed lines correspond 
to approximations obtained from the linear stability analysis. 
The arrows indicate the time flow. 



clearly breaks down very quickly, since the linear stability 
analysis shows an unstable spiral, which continues to only 
grow about T+. 

At an even larger 5 value, we see in Fig. 4(c) that the 
numerical solution immediately diverges from the fixed 
point and goes to an oscillatory state, the amplitude of 
which is limited by the nonlinear damping term in Eq. 
(2). 

The behavior for other values of a is analogous. There- 
fore, we have identified within the DAO model (2) a 
whole region of a and S values where we obtain oscil- 
lations of the SST anomaly between a hot (T + > 0) and 
a cold (T_ < 0) fixed point. This feature is what makes 
the DAO a model for oscillations between hot El Nino 
and cold La Nina conditions. 



E. Comparing periodicities 

As identified in the last section, the DAO model is 
capable of producing limit cycle oscillations for a wide 
region of (a, S) values. We next compute the period of 
these oscillations and compare them to those observed 
for typical ENSO cycles. 
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As shown in Figure 2, the Nino 3.4 region extends 
from 120° West to 170° West, with the midpoint at 145° 
West. If we assume that the Pacific western boundary is 
at w 120° East, this gives an angular separation of 95° 
of longitude for the waves to travel. This corresponds to 
a distance of 95(27r/360)rEarth = 10.6 x 10 6 m, where the 
radius of the Earth rEarth = 6.37 x 10 6 m. Using the val- 
ues of lims^ 1 for a Kelvin wave and 0.47ms _1 for the 
primary Rossby mode as obtained in Appendix A thus 
gives a delay of 262 days for the Rossby propagation to 
the western boundary and a further 87 days for the re- 
turn of the Kelvin waves, a total delay of A = 349 days. 

This information now allows us to compute the time 
scaling parameter, however since k = 5/A, these are dif- 
ferent for each value of S used in the DAO model. By the 
use of a Fourier Transform applied over 600 years, we 
find e.g. for a = 0.7, 6 = 3 a period r' = 11.1 which cor- 
responds to a true DAO period of r = r'/fc = 3.5 years. 
In Fig. 6 we show the results of many runs with varying 
a = 0.56, 0.57, . . . , 1 and 5 = 1, 1.1, . . . , 4.9. Comparing 




FIG. 6: Grey scale plot of DAO periods for various a and <5 
values. Lines of constant periods are indicated, those corre- 
sponding to periods 3.3 and 3.7 are shown as dashed lines. 
The white region corresponds to the stable region as present 
in Fig. 3, with the thick black line representing the first neu- 
tral curve. 

these periods to real ENSO events, which have a period 
of 3.3 to 3.7 years, we see that there is a large region of 
(a, S) values where there is a good agreement. 

F. Discussion of the basic DAO model 

In the preceding sections, we have seen that the DAO 
is a useful tool in studying the ENSO cycle. The main 
beauty of the DAO (1) is its simplicity. Although the 
full equation cannot be solved explicitly, the linear stabil- 
ity analysis and the numerical solutions both give highly 
useful insights. In particular, the numerical solutions ex- 
hibit the true global dynamics of the DAO. There are 



limit cycle oscillations, with a period that matches the 
real world, and hence the DAO is a good model for these 
features of ENSO. Nevertheless, the DAO model assumes 
a localized coupled problem, which is obviously an over- 
simplification of the situation in the Pacific. All influ- 
ences of the atmosphere, the earth's rotation, and the 
ocean currents and dynamics are simply modeled by a 
non-linear, effective damping term and a somewhat ar- 
bitrary delay. The delay term ignores, e.g. East bound- 
ary reflections that could possibly be considered. The 
model further assumes that time scales must be related to 
propagation times of equatorially trapped oceanic waves 
in a closed basin. For the atmosphere, which has short 
time scales, the DAO simply assumes an instantaneous 
response. 



IV. ANNUAL FORCING 

In this and the following sections, we shall extend the 
basic DAO discussed above by including external influ- 
ences which allow us to model the ENSO dynamics bet- 
ter. In order to distinguish this extended DAO from its 
original, we shall denote the new model as ENSO-DAO 
from now on. 



A. Modifying the DAO to include the annual cycle 

As mentioned above, an El Nino event typically ap- 
pears in late December, early January. Thus the annual 
cycle is clearly playing an important role in determining 
the onset of El Nino [18, 19]. Let us thus add to the 
DAO (1) another term that models the periodic heat- 
ing and cooling of the equator during the annual cycle. 
The average SST data, as shown on Fig. 1, can be used 
to construct a Mathematica interpolating function, Y(t), 
representing the continuous changes of the average SST 
during the year [27]. We therefore include its numerical 
time derivative, representing the expected change in SST 
due to annual forcing, to give a new DAO equation 

^ = kT- bT 3 - AT(t - A) + ^Y(t) (9) 

(XL CLl 

This equation now depends critically on the (a, S) values 
chosen in the dimcnsionlcss model, since computation of 
the scaling constants k and b is required. We have seen 
earlier that k = 5/A, and b is given by the expression 
b = k(T'/T) 2 . In order to estimate b we shall take T'/T 
as the ratio of their respective maxima/minima. Figure 
1 shows the maximum anomaly is on average 2°C, whilst 
the maximum value of T" can be calculated for any values 
(a, S). 

Choosing a = 0.7, 5 = 3 as the standard ENSO-DAO 
model for the next few sections allows us to calculate 
k = 3.14/years, b = 1.09 °C~ 2 /years, and hence discuss 
the behavior of this specific model under annual forcing. 
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B. Solving the annually forced DAO 

Fixed points of (9) do not exist; they are now fixed 
cycles fluctuating about 0, ±y/(k — A)/b due to the pe- 
riodic nature of the forcing term. Nevertheless, taking 
t = to be the start of the year, the expected SST will 
be at its coldest and dY/dt = 0, so we can take the initial 
condition to be T(0) = y/ik — -4)/& + 0.15 (this is similar 
to before but now includes the scaling constants). 

Let us now consider numerical solutions of Eq. (9). In 
Fig. 7, we show the solution of T vs t for parameters as 
above. The annual cycle as in Fig. 1 is also shown. We 
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FIG. 7: Numerical solution of the annually forced DAO (9) 
with 5 = 3 and a = 0.7 showing the annual cycle Y(t) (dashed 
black), model solution T (grey) and SST anomaly T — Y(t) 
(solid black). 

note that with the inclusion of annual forcing, T is no 
longer the SST anomaly, but rather, the anomaly is the 
difference between T and the annual forcing, Y(f). This 
true SST anomaly is also shown in Fig. 7, and should be 
compared to Fig. 1. 

From Fig. 7 it is clear that the addition of an annual 
forcing term has changed the dynamics in several ways. 
Firstly, the peak amplitude of the oscillations has been 
increased. The peaks themselves are less regular because 
of the quickly changing forcing term provided by the sea- 
sons. However, the paramount effect of the annual cycle 
is in determining the onset and period of an El Nino 
event. Hot events only occur in-phase with the annual 
cycle, and hence the ENSO-DAO now has an exact inte- 
ger period, e.g. 3 years for our model. The ENSO-DAO 
gives high temperatures from July through to March of 
the next year — which is indeed as observed in the Pacific 
— with the maximum amplitude usually in December. 



V. GLOBAL WARMING 

Global warming is another effect that may influence 
the periodicity and amplitude of ENSO events. For the 
purpose of this paper it will be modeled via a constant 



heating on the Pacific DAO system. Working with the 
dimcnsionless DAO for simplicity, this is achieved by 
adding a constant (3 to the DAO equation (2) giving: 



-r= T 
dt 



aT(t-S)+ (3 



(10) 



On its own, this term would give ^rr = P, solved by 
T = (3t+c, a linear warming with time. We shall consider 
(3 > 0, a constant heating, although (3 < cooling can 
be treated completely analogously. 

Current climate change models predict the rise in av- 
erage global temperature over the next 100 years will be 
below 5°C [20], hence we chose this scenario to be used 
in our ENSO-DAO. Using the standard model of Section 
IV, this gives 



b 5°C 

warming - y fc IQQyearS 



/^global 

as the value of (3 in Eq. (10). 



0.009 



(11) 



A. Fixed points, linear stability analysis and 
consequences 

The fixed point solutions are now determined by the 
equation: 



(l-a)T(«-(T(«) 3 + /? = 







(12) 



For small (3, there is hardly any change from the original 
T fixed points, and hence the system will behave as the 
basic DAO. As (3 is increased, the unstable inner fixed 
point T (/3) dc creases from zero, whereas both the outer 

solutions T± increase. In cases when (3 is large, and for 
large enough a, only one fixed point remains on the real 
axis (the other two are now complex). 

Instead of considering solutions starting near T + = 
\f\ — a, we now start near . The neutral curves can 
be determined via the same methods used to obtain Eq. 
(7), but with the fixed point f'{ ,3) now depending on [3: 



2nn =F arccos 


"l-3(Tf >) 2 " 


a 


\l 




i-3 (if ) 2 " 


2 



n e NU {0} (13) 



The first neutral curve for a /3 g i bai warming DAO system, 
compared to no heating, is plotted in Fig. 8. 

Consider a basic DAO system (above the solid black 
neutral curve in Figure 3) that has an unstable fixed 
point ff -0 ^ with limit cycle oscillations. By introduc- 
ing (3 > 0, the neutral curve is shifted to the right, and 
hence the additional (hashed) stable region in Fig. 8 is 
created, changing the DAO dynamics. As (3 increases, 
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FIG. 8: The first neutral curve for the original DAO (solid 
black, as in Fig. 3), forced DAO (dashed black) and symmetri- 
cally coupled DAO (dotted black) . Also shown is the stability 
region affected by global warming (hashed), and by coupling 
(dotted). Higher mode neutral curves are not shown. 

the model fluctuates for longer about the unstable T_j^ 
before nipping out to oscillate around both fixed points. 
Eventually, the dashed neutral curve in Fig. 8 moves to 
the right of (a, 5) and 7+ becomes attracting. This sta- 
bility bifurcation occurs at /?bif, which can be found for 
any given (a, 6) by solving Eq. (13) for (3. 

If enough heating were applied to the DAO, the ENSO 
cycle would be terminated, leaving an ocean at a con- 
stant temperature and with no El Nino events. However, 
for slow heating, such as (5 — 0.009, there is very little 
change to the time dependence of the DAO with respect 
to SST anomaly amplitudes or their periodicity except 
for a small subset of solutions. This would suggest that 
global warming will not affect ENSO, and hence El Nino 
events will still continue to occur. 



VI. STOCHASTIC EFFECTS 

Chaos inherent in real world weather systems plays an 
important part in the irregularity of El Nino events. Pro- 
cesses such as storms, cyclones and ocean currents which 
traverse the Pacific play no part in our ENSO model, but 
will influence the system for their duration in the Nino 
3.4 region. We therefore consider the addition of a gen- 
eral stochastic term [21] R(t) to the dimensional DAO 
equation (1) giving: 

AT 

— = kT - bT 3 - AT(t - A) + R(t) (14) 

R(t) is a step function that takes fixed, normally dis- 
tributed (mean 0, standard deviation A), random val- 
ues for each monthly time-step. A (non-ENSO) event 
that lasts for longer than one month is unphysical, and 
shorter-lived events will not have time to cause any sig- 
nificant differences (their randomness is averaged out). 



Another way of viewing R(t) is as a j3 heating/cooling 
that only lasts for one month. This will mean that the 
neutral curves given by Eq. (13) will shift right/left each 
month, due to the fixed points Tj_ changing position. 

If A is small, then there is hardly any change to the ba- 
sic DAO. However, for large enough A, the fixed points 
can change their stabilities each month, inducing a great 
deal of random behavior in the DAO. For example, the 
system may be stable and tending towards of that 
month, and then R(t) changes to produce an unstable 
system that wants to oscillate. This creates many differ- 
ent El Nino events with differing amplitudes and dura- 
tions, as well as the possibility of times where the system 
stays near zero. 



VII. THE FULL ENSO-DAO 

We now combine of all of the effects considered above, 
namely the annual cycle, global warming and stochastic 
effects. Therefore the DAO equation (1) becomes 

AT 

— = kT- bT 3 - AT(t-A)+ Y'{t)+ B+ R(t) (15) 

Using the standard model of Section IV, with B = 
5°C/100years as the fully dimensional global warming 
factor and R(t) with A = 5, the DAO produces very real- 
istic pictures, such as Fig. 9, that look remarkably sim- 
ilar to Fig. 1. Most obvious are the small fluctuations 

4 [ ' Combined DAO Solution ' 1 ' I 
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FIG. 9: Numerical solution of the combined DAO (15) with 
5 = 3 and a — 0.7 showing the annual cycle Y(t) (dashed 
black, right axis), model solution T (grey, left axis) and SST 
anomaly T — Y(t) (solid black, left axis). 

associated with the stochastic effects. In addition, the 
regularity of the original DAO oscillations is broken and 
sometimes a complete ENSO event has vanished. E.g. 
the absence of the SST peaks at about t <~ 20 in Fig. 9. 
This happens when the onset of the ENSO cycle event 
takes place at a time when the random forcing happens 
to change the SST into the opposite direction. We em- 
phasize that a variety of solution curves similar to Fig. 
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9 can be generated by taking other (a, 6) combinations, 
initial SST values and stochastic forcing. 



VIII. COUPLED DAOS 

The basic DAO model (1) considers a single, but repre- 
sentative region with strong atmospheric-ocean coupling 
in the Pacific. Clearly, we might also want to consider 
situations in which other regions of the Pacific are in- 
corporated. These could be either other regions along 
the equator, where we would expect differing coupling 
strengths and delay times. Or they could be regions to 
the North and South of the central DAO region in which 
no delay and only weak coupling exists, such that these 
additional regions act as temperature sinks. Possible ar- 
rangements are shown in Fig. 10. In this chapter, we 
will consider one such prototype situation, namely two 
coupled regions along the equator. 



N 
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The fixed point solutions can be found as before by 
taking = and substituting (16a) into (16b). This 

yields a 9th order polynomial in with a mixture of 
real and complex roots which vary as a function of 7. 

The solutions 7^ = ±\A — a + 7 are the same in both 
regions and are always real since < a < 1, and hence 
are of most interest. 

Conducting a linear stability analysis about the fixed 
points (T"i+ , ?2+) as in Section IIIC, we find neutral 
curves given by 

2n7r T arccos (^^27) r , 

5= i V g - , neNU{0} (17) 

yW 2 - (3a - 2 - 2 7 ) 2 

The first neutral curve is shown in Figure 8 for an exam- 
ple with 7 = 0.2. 

We see that similar to Section V, the neutral curves 
move right with increasing 7, such that the fixed point 
can bifurcate from unstable to stable. So, if there is 
strong coupling between regions in the Pacific (large 
7), then independent oscillations cannot occur and the 
system settles down to constant temperature anomalies. 
However, for sufficiently weak coupling, ENSO-style os- 
cillations are still permitted. 

As before, we have also solved the Eq. (16) numerically. 
Depending on the initial SST anomalies, we find that we 
can drive the regions into (i) stable, non-oscillatory solu- 
tions, (ii) in-phase oscillations and even (iii) oscillations 
which are out-of-phase by it. For brevity, we do not show 
the results here. 



FIG. 10: Schematic of the Pacific showing the original ENSO- 
D AO region (grey) , new coupled regions 1 and 2 and potential 
North-South regions (hashed). The equator is indicated as a 
horizontal dashed line. 



A. Two coupled DAOs along the equator 

Let us for simplicity assume that coupling a and de- 
lay S in the two regions are the same. This is a valid 
assumption as long as the two adjacent regions are both 
in the strong-coupling region such that the distance to 
the Western boundary is approximately the same, with 
the same losses and reflection properties. If one region 
is warmer than the other, then the flow of heat energy 
across their common boundary will act to heat one and 
cool the other. This can be modeled by incorporating 
a "fTi + i term in each DAO with 7 > and Tj+i the 
anomaly in the adjacent DAO region. Then we have a 
system of 2 coupled DAOs such that 

rIT 

—j^- = Ti — T x 3 - oTi (t - S) + 7 T 2 (16a) 
cIT 

— - = T 2 - T 2 3 - aT 2 (t -S) + 7T1 (16b) 



B. Two regions with different ENSO coupling 

Let us now consider briefly the situation in which a, or 
Si arc different in the different regions. If T 2 represents 
the anomaly in region 2 of Fig. 10 and Ti the anomaly 
in region 1, then the wave transit time from Tl to the 
western boundary would be greater than that from T 2 . 
Similarly, the strength of the returning signal to T\ would 
be less than that to T 2 , as more information would be lost 
due to the greater distance across which information is 
transported. However, the waves propagating from T\ 
are not hidden from the problem whilst they are propa- 
gating through T 2 , suggesting that S could be the same 
for both regions, namely the propagation time from T 2 to 
the western boundary and back. Similarly, the strength 
of the local coupling is assumed to be less in T\ , suggest- 
ing that even though the strength of the returning signal 
is less, relative to the local coupling its affect would be the 
same as in T 2 . Figure 11 shows how the coupled model 
(16) is affected by different values of a. In the Figure we 
see that T\ shows some irregular behavior with varying 
amplitudes. This behavior is more like that of El Nino 
than the simple DAO. 

Other possible coupled DAO scenarios could include 3 
coupled regions, either along the equator or aligned along 
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FIG. 11: Numerical solutions to the coupled DAO (16), with 
ai = 0.5 in region 1, Q2 = 0.75 in region 2, 8 = 4 in both 
regions and 7 = 0.1. 



a south-north axis. Our preliminary results indicate that 
such models gives rise to irregular ENSO cycles as above. 
However, the models have to be treated with proper care, 
since some solutions become unphysical, e.g. indicating 
uninterrupted heating or cooling of the ocean. 



IX. CONCLUSIONS 

The DAO is given by a simple one-dimensional differ- 
ential equation, and may nowadays be studied by readily 
available desk- or laptop computing facilities. Yet it cap- 
tures much of the observed richness and variance of the 
ENSO phenomenon. We are therefore convinced that it 
presents an ideal tool to introduce students to the topic. 
Let us briefly recapitulate what issues we have been able 
to study in the present paper using the DAO. We first 
outline the linear stability analysis and the numerical so- 
lution of the DAO that were the content of the first paper 
on the DAO [8]. This then serves to introduce the con- 
cept of shallow-water Kelvin and Rossby waves, as well as 
to establish our notation. Using the corresponding wave 
speeds and the geography of the Pacific, we can already 
show that for a large set of parameters, the observed 
ENSO periodicity can be modeled by the DAO. 

Next, the inclusion of the measured annual cycle 
demonstrates that the observed ENSO-DAO periodicity 
may in fact be heavily influenced by the Earth's period 
around the Sun. The occurrence of an ENSO event is 
very much subject to the right season and thus the orig- 
inal DAO period need not be identical to the observed 
ENSO-DAO period. We stress that the modeling of the 
annual cycle within the DAO is quite easy and should 
not be beyond any student who has mastered the origi- 
nal DAO. 

We then extended the DAO to include a constant heat- 
ing — global warming — and found that current pre- 
dictions of 5°C rise over 100 years are unlikely to in- 



fluence the ENSO-DAO dynamics. We also do not find 
any large increase in either occurrence or magnitude of 
ENSO events. This seems to be at variance with current 
predictions based on CGCMs [22, 23] and might just in- 
dicate the limitations of the ENSO-DAO model for true, 
quantitative predictions [28]. 

Putting it all together, we then included a certain 
level of stochastic forcing alongside the annual cycle, and 
found that the ENSO-DAO dynamics now become very 
much like the observed ENSO dynamics, showing a large 
variation in oscillation periods and amplitudes. At least 
on a qualitative level — the range of the SST oscillations 
is in fact even quantitatively similar — which shows the 
strength of the ENSO-DAO. 

Finally, we briefly discussed some extensions of the 
ENSO-DAO model, such as coupled regions. The ex- 
ample of the 2-couplcd-regions model shows that the 
extended DAOs share some irregularities with the ob- 
served ENSO phenomenon already, even in the absence 
of annual or stochastic forcing. Further investigations of 
ENSO-DAO might include an analysis of the statistical 
properties of the measured ENSO cycles as in Fig. 9 and 
comparisons with observed data, a more thorough anal- 
ysis of the possible changes in periodicity and strengths 
of ENSO-DAO oscillations due to human influences, a 
mathematical analysis of the DAO for integer periods or 
even an application to the Atlantic. 

In conclusion, we see that the DAO — while being a 
research tool in its own right [9] — can be investigated 
by modest analytical and computational resources. We 
think that it will be an ideal topic for all students, teach- 
ers and lecturers interested in the exciting dynamics of 
the El Nino phenomenon and hope that the present paper 
will help them in this quest. 
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APPENDIX A: SHALLOW WATER EQUATIONS 
AND WAVES 

1. Two-fluid model and shallow water equations 

Oceans are warmer near the surface than they are in 
the deep, and due to thermal expansion this warm layer 
will have a lower density. Thus a simple model for the 
ocean structure is a two-fluid model with density p\ in 
the warm surface layer and density pi > p\ below in the 
cold abyss. The dividing zone is called the thermocline. 
At its widest, the Pacific is 17, 700km across from the 
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Malay Peninsula to Panama whilst only having an aver- 
age depth of 4, 280m. So the Pacific is roughly 3 orders of 
magnitude wider than it is deep, and so we are justified in 
rewriting the original Navier-Stokes equations using the 
shallow-water approximation [24, 25]. These are given as 



drj ^ ( du dv 

dt \ dx dy 

du du du .dri 

at ox dy ox 

dv dv dv dn 

77T + u 77~ + + f u + -9 7T 

dt dx dy dy 



(Ala) 
(Alb) 
(Ale) 



in the absence of external forces. Here (u, v, w) is the dis- 
placement velocity in the Cartesian frame (x,y,z) and 
r](x,y,t) is a small perturbation to the constant mean 
thermocline depth H. Furthermore, / is the Coriolis pa- 
rameter and g' — P2 ^ Pl g is the effective gravitational ac- 
celeration. Equations (Al) permit two types of thermo- 
cline wave solutions. The first are inertia-gravity Kelvin 
waves. The second type are planetary Rossby waves. 



2. Kelvin waves 

Equatorially trapped Kelvin waves have no meridional 
(y) velocity fluctuations such that the shallow water 
equations (Al) can be simplified to 



d 2 u 'n^ 2u o 
dt 2 ^ dx 2 



f 



du 
~dt 



9'H 



d 2 u 
dxdy 



(A2a) 



(A2b) 



Equation (A2a) implies that solutions have to be of the 
form u = E{y)F(x ± ct) with wave speed c = y/g'H. 
In addition, Eq. (A2b) prohibits westward propagation 
{F{x + ct)) since those solutions become unbounded at 
large y. Therefore only eastward, equatorially trapped 
(v = 0), Kelvin waves are possible. These non-dispersive 
thermocline gravity waves propagate eastward with speed 
c w 1.4ms -1 , when assuming a height H = 100m and 
(P2 - Pi)/ Pi = 0.002 [11] for the two-fluid model [29]. 



3. Rossby waves 

Rossby waves owe their existence to the variation of the 
Coriolis parameter with latitude. Close to the equator 
the Coriolis parameter is given by / = 2f2y/rEarth where 
f Earth is the radius of the Earth and the Earth's angular 
velocity. If we consider a linearized version of the shallow- 
water equations and seek solutions for latitudinal waves 
of the form v = V(y) expi(fcx — u>t), with u> > and the 



sign of k determining the direction of zonal [x) phase 
propagation, we obtain their wave equation as 



d 2 v p 2 

dy 2 + g'H 



(Y 2 -y 2 )V = 



with 



d 2 {g'H 



(A3) 



(A4) 



and (3 = 2fi/rEarth- From (A3) we see that solutions are 
wavelike within an equatorial waveguide of width 2Y, and 
decay exponentially for latitudes greater than \Y\. 

If we now assume that the latitudinal dependance of V 
is of plane- wave type ~ exp(int/), then for these equato- 
rially trapped waves we have n € N as the only allowed 
standing wave solutions. A dispersion relation can be 
extracted from (A3) and at low frequencies 



-(3k 



k 2 + 2n+l 



n € N 



(A5) 



with wavelength A = \J c/ (3. Since k must therefore be 
negative, all Rossby waves have westward phase propa- 
gation. The slow, short, dispersive waves have eastward 
group velocities and the fast, long, non-dispersive waves 
have westward group velocities of c/(2n + 1). Only the 
long waves are important in the oceanic adjustment to 
changes in surface winds — a process that is crucial in 
ocean-atmosphere coupled systems — so only westward 
propagating Rossby waves will usually be considered. 

The fastest Rossby wave (n = 1) travels at one third 
the speed of a Kelvin wave, e.g. w 0.47ms _1 in the Pa- 
cific. These are the non-dispersive, equatorially trapped, 
westward propagating Rossby waves considered in the 
DAO model in Section III. 



APPENDIX B: NUMERICAL ROUTINES 

We include an abbreviated, but fully functional Mafch- 
ematica 5.2 code for the solution of the annually forced 
DAO with global warming and monthly random weather 
as studied in Section VII. Physical units used in the code 
are given in years for time and Kelvin for temperature 
throughout. 

«NDelayDSolve.m 

(* Define the variables to be used and calculate 
the required scalings *) 

year=365.24; delay=349/year ; alpha=0.7; delta=3; 
beta=0.05; k=delta/delay ; time=25; 

(* Solve the basic dimensionless DAO, determine its 
maxima and therefore the value of b, the 
temperature scaling *) 

Tf p=Sqrt [1-alpha] ; 
orig= NDelayDSolve [ 
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-CT'[t] == T[t]-T[t]-3-alpha*T[t-delta]}, 
{T->(Tfp+0.05&)}, 
{t,0,25*delta}, 
MaxRecursion->50 
]; 

values=Table [T [t] / . orig , {t , , 25*delta ,0.1}]; 

Tmax=Max [values] ; 

b=k*(Tmax/2)-2; 

(* Read the real data from file, and construct 
annual cycle function *) 

data=ReadLlst ["ave_sea_temp . txt "] ; 
yearly = Interpolation [data-27 . 1 , 

PeriodicInterpolation->True] ; 
Y[t_] :=yearly [ 12*t +1 ] 

O Create a Gaussian Distributed random variable 
(mean zero, standard deviation lambda) which 
changes every month, representing the effects 
of stochastic processes on the DAO *) 

lambda=5 ; 

Needs ["Statistics 'ContinuousDistr ibut ions ' "] ; 
stochastic=Table [Random [ 

NormalDistribution [0, lambda] ] , {time*12+l} 

]; 

R[t_] := stochastic [ [IntegerPart [12*t] +1] ] ; 
Plot [R[t] ,{t,0,time}] 

(* Solve the Combined DAO Equation *) 

Tf = Extract [T /. 

Solve [k*T - b*T~3 - k*alpha*T + beta == 0, T] , 
3]; 



combo = 
NDelayDSolve [ 
{(T')[t] == k*T[t] - b*T[t]~3 
- k*alpha*T[t - delta/k] 
+ beta + R[t] + (Y') [t]}, 
{T->((Tf + 0.15 &)) 
}, 

{t, 0, time}, 

MaxRecursion -> 100, MaxSteps -> Infinity 

]; 

(* Plot the Combined DAO Solution, the Annual Cycle, 
and the SST-DA0 Anomaly *) 

Plot[ 

{Evaluate [T [t] / . combo] , Y [t] , Evaluate [T [t] / . combo] 
-Y [t]},{t,0,time}, 
AxesLabel->{"t/years" , "T/K"}, 
PlotStyle->{{RGBColor [0,1,0]}, 
{RGBColor [0,0,0]}, 
{RGBColor [1,0,0]}} 

]; 

(* Calculate the period of the anomaly 

oscillations using Fourier Transform method *) 

data=Abs [Fourier [Table [ (T [t] / . combo) -Y [t] , 

-ft, 0, time, 0.1}]]] ; 
Do [If [data [ [k] ] ==Max [data] , ans=Return [k] ] , 

{k,l,time,l}] ; 
period=N [t ime/(ans-l)] 
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